////////////////////////
// L为下三角矩阵，Ly=b
// n为矩阵的秩
// 解得的y储存在b中
////////////////////////
int Lower_Triangular_Matrix(double** L, double* b, int n);
////////////////////////
// L为单位下三角矩阵，Ly=b
// n为矩阵的秩
// 解得的y储存在b中
////////////////////////
int Unit_Lower_Triangular_Matrix(double** L, double* b, int n);
////////////////////////
// U为上三角矩阵，Ux=y
// n为矩阵的秩
// 解得的x储存在y中
////////////////////////
int Upper_Triangular_Matrix(double** U, double* y, int n);
////////////////////
// A为n×n矩阵，求A的三角分解
// A=LU，其中L为单位下三角阵，U为上三角阵
// 结果储存在A的内存中，A的上三角部分为U，其余部分为L除主对角线外的元素
////////////////////
int Gaussian_Elimination(double** A, int n);
//////////////////////////////
// 列主元高斯消去法Ax=b
// A为n×n矩阵，求A的三角分解
// PA=LU，其中L为单位下三角阵，U为上三角阵
// 结果储存在A的内存中，A的上三角部分为U，其余部分为L除主对角线外的元素
//////////////////////////////
int Column_Principle_Gaussian_Elimination(double** A, double* b, int n);




///////////////////////////////////////////////
#include <iostream>
#include <math.h>



////////////////////////
// L为下三角矩阵，Ly=b
// n为矩阵的秩
// 解得的y储存在b中
////////////////////////
int Lower_Triangular_Matrix(double** L, double* b, int n)
{
  for(int j = 0; j < n-1 ; j++)
    {
      b[j] = b[j]/L[j][j];
      for(int i = j + 1; i < n ; i++)
	{
	  b[i] = b[i]-b[j]*L[i][j];
	}
    }
  b[n-1] = b[n-1]/L[n-1][n-1];
  return 0;
};

////////////////////////
// L为单位下三角矩阵，Ly=b
// n为矩阵的秩
// 解得的y储存在b中
////////////////////////
int Unit_Lower_Triangular_Matrix(double** L, double* b, int n)
{
  for(int j = 0; j < n-1 ; j++)
    {
      for(int i = j + 1; i < n ; i++)
	{
	  b[i] = b[i]-b[j]*L[i][j];
	}
    }
  return 0;
};



////////////////////////
// U为上三角矩阵，Ux=y
// n为矩阵的秩
// 解得的x储存在y中
////////////////////////
int Upper_Triangular_Matrix(double** U, double* y, int n)
{
  for(int j = n-1; j > 0 ; j--)
    {
      y[j] = y[j]/U[j][j];
      for(int i = 0; i < j ; i++)
	{
	  y[i] = y[i]-y[j]*U[i][j];
	}
    }
  y[0] = y[0]/U[0][0];
  return 0;
};



////////////////////
// A为n×n矩阵，求A的三角分解
// A=LU，其中L为单位下三角阵，U为上三角阵
// 结果储存在A的内存中，A的上三角部分为U，其余部分为L除主对角线外的元素
////////////////////
int Gaussian_Elimination(double** A, int n)
{
  for(int k = 0; k < n - 1; k++)
    {
      for(int i = k + 1; i < n; i++)
	{
	  A[i][k] = A[i][k]/A[k][k];
	}
      for(int i = k + 1; i < n; i++)
	{
	  for(int j = k + 1; j < n; j++)
	    {
	      A[i][j] = A[i][j] - A[i][k]*A[k][j];
	    }
	}
    }
  return 0;
};

//////////////////////////////
// 列主元高斯消去法息
// A为n×n矩阵，求A的三角分解
// PA=LU，其中L为单位下三角阵，U为上三角阵
// 结果储存在A的内存中，A的上三角部分为U，其余部分为L除主对角线外的元素
//////////////////////////////
int Column_Principle_Gaussian_Elimination(double** A, double* b, int n)
{
  for(int k = 0; k < n - 1; k++)
    {
      int p = k;   //p为第k列中最大的元素行数
      for(int i = k + 1; i < n; i++)
	{
	  if(fabs(A[i][k]) > fabs(A[p][k]))
	    p = i;
	}
      if(p != k)   //对A进行初等变换，把k、p行互换
	{
	  for(int i = 0; i < n; i++)
	    {
	      double tmp = A[p][i];
	      A[p][i] = A[k][i];
	      A[k][i] = tmp;
	    }
	  double tmp = b[p];
	  b[p] = b[k];
	  b[k] = tmp;
	}
      
      if(A[k][k] != 0)
	{
	  for(int i = k + 1; i < n; i++)
	    {
	      A[i][k] = A[i][k]/A[k][k];
	    }
	  for(int i = k + 1; i < n; i++)
	    {
	      for(int j = k + 1; j < n; j++)
		{
		  A[i][j] = A[i][j] - A[i][k]*A[k][j];
		}
	    }
	}
      else
	{
	  std::cout << "Error! A is a singular matrix." << std::endl;
	  return -1;
	} 
    }

  return 0;
};




